///\todo upgrade

//////////////////////////////////////////////////////////////////////////////////
//                                                                              //
//  This file is part of the buola project (https://code.google.com/p/buola/).  //
//                                                                              //
//  Copyright(c) 2007-2012 Xavi Gratal                                          //
//  gratal AT gmail DOT com                                                     //
//                                                                              //
//  Buola is free software: you can redistribute it and/or modify               //
//  it under the terms of the GNU General Public License as published by        //
//  the Free Software Foundation, either version 3 of the License, or           //
//  (at your option) any later version.                                         //
//                                                                              //
//  Buola is distributed in the hope that it will be useful,                    //
//  but WITHOUT ANY WARRANTY; without even the implied warranty of              //
//  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the               //
//  GNU General Public License for more details.                                //
//                                                                              //
//  You should have received a copy of the GNU General Public License           //
//  along with buola.  If not, see <http://www.gnu.org/licenses/>.              //
//                                                                              //
//////////////////////////////////////////////////////////////////////////////////

#ifndef _BUOLA_MAT_MARKOV_H_
#define _BUOLA_MAT_MARKOV_H_

#include <buola/mat/cmat.h>
#include <buola/mat/initializers.h>

namespace buola { namespace mat {

CVec_i generate_mm(size_t pLen,const CMat_d &pT)
{
    assert(pLen>0);
    assert(pT.Rows()==pT.Cols());
    
    size_t lN=pT.Rows();
    
    mat::CVec_i lX(pLen);
    lX(0)=0; //we assume we are originally in state 0
    
    for(int i=1;i<pLen;i++)
    {
        const auto &lPRow=pT(lX(i-1),nAll);
        
        double lCQX=max(lPRow);
        while(true)
        {
            lX(i)=random::get<int>(lN-1);
            double lU=random::get<double>();
            if(lU<lPRow(lX(i))/lCQX) //only accept it if u<p_x/(c*q_x)
                break;
        }
    }
    
    return std::move(lX);
}

CMat_d estimate_mm(const CVec_i &pStates,size_t pN)
{
    //compute the frequency for each transition and each state
    mat::CMat_d lF=mat::zeros(pN,pN);

    for(int i=1;i<pStates.Rows();i++)
        lF(pStates(i-1),pStates(i))++;
    
    lF/=extend(sum(rows(lF)));

    return std::move(lF);
}

std::pair<mat::CVec_i,mat::CVec_i> generate_hmm(size_t pLen,const mat::CMat_d &pA,const mat::CMat_d &pB,
                                                  const mat::CRow_d &pPi)
{
    mat::CVec_i lX(pLen),lY(pLen);
    
    mat::CRow_d lARow=pPi;

    for(int i=0;i<pLen;i++)
    {
        //generate X sample according to distribution
        lX(i)=sample(lARow);
        //select observation probabilities given matching X
        lY(i)=sample(pB(lX(i),nAll));
        
        lARow=pA(lX(i),nAll);
    }
 
    return {lX,lY};
}
#if 0
mat::CVec_i viterbi_hmm(const mat::CMat_d &pA,const mat::CMat_d &pB,const mat::CRow_d &pPi,const mat::CVec_i &pY)
{
    int lN=pY.size();
    int lM=pPi.size();
    
    mat::CMat_d lD=mat::zeros(lN,lM);
    mat::CMat_d lU=mat::zeros(lN,lM);
    
    lD(0,nAll)=log(pPi);
    mat::CMat_d lLogA=log(pA);
    mat::CMat_d lLogB=log(pB);

    for(int i=1;i<lN;i++)
    {
        for(int j=0;j<lM;j++)
        {
            //check from which previous state it is 'cheaper'
            //to get to state j. Keep the index in u
            lD(i,j)=max(lD(i-1,nAll)+lLogA(nAll,j).T(),lU(i,j));
            //take observation probability into account
            lD(i,j)+=lLogB(j,pY(i));
        }
    }
    
    mat::CVec_i lX=mat::zeros(lN,1);
    lX(lN-1)=indmax(lD(lN-1,nAll));
    for(int i=lN-2;i>=0;i--)
        lX(i)=lU(i+1,lX(i+1));
    
    return lX;
}
#endif

/*namespace mat*/ } /*namespace buola*/ }

#endif
